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Abstract 

We consider several systems of nonlinear hyperbolic conservation laws describing the dynam- 
ics of nonlinear waves in presence of phase transition phenomena. These models admit under- 
compressive shock waves which arc not uniquely determined by a standard entropy criterion but 
must be characterized by a kinetic relation. Building on earlier work by LeFloch and collabora- 
tors, we investigate the numerical approximation of these models by high-order finite difference 
schemes, and uncover several new features of the kinetic function associated with with physi- 
cally motivated second and third-order regularization terms, especially viscosity and capillarity 
terms. 

On one hand, the role of the equivalent equation associated with a finite difference scheme is 
discussed. We conjecture here and demonstrate numerically that the (numerical) kinetic function 
associated with a scheme approaches the (analytic) kinetic function associated with the given 
model — especially since its equivalent equation approaches the regularized model at a higher 
order. On the other hand, we demonstrate numerically that a kinetic function can be associated 
with the thin liquid film model and the generalized Camassa-Holm model. Finally, we investigate 
to what extent a kinetic function can be associated with the equations of van der Waals fluids, 
whose flux-function admits two inflection points. 

Key words: hyperbolic equation, conservation law, shock wave, kinetic relation, viscosity, capillarity, 
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1. Introduction 

1.1. Background 

In this paper wc study the numerical approximation of several first-order nonlinear 
hyperbolic systems of conservation laws, and we consider discontinuous solutions gener- 
ated by supplementing the hyperbolic equations with higher-order, physically motivated, 
vanishing rcgularization terms. Specifically, wc consider complex fluid flows when phys- 
ical features such as viscosity and capillarity effects can not be neglected even at the 
hyperbolic level of modeling, and need to be taken into account. This gives rise to many- 
shock wave theories associated with any given nonlinear hyperbolic system: depending 
on the underlying small-scale physics of the problem under consideration, one need a 
different selection of "entropy solutions". 

It is well known that solutions of nonlinear hyperbolic equations become discontinuous 
in finite time, and may, therefore, exhibit shock waves. Classical shock waves satisfy 
standard entropy criteria (due to Lax, Oleinik, Wendroff, Liu, Dafermos, etc); they are 
compressible and stable under perturbation and approximation. 

On the other hand, discontinuous solutions of hyperbolic problems may also exhibit 
non-classical, under-compressive shock waves - also referred to as subsonic phase bound- 
aries in the context of phase transition theory. To uniquely characterize under-compressive 
shocks one need to impose a jump condition that is not implied by the given set of con- 
servation laws and is called a kinetic relation. The selection of physically meaningful 
shock waves of a first-order hyperbolic system is determined by traveling waves associ- 
ated with an augmented system that includes viscosity and capillarity effects. That is, 
one searches for scale-invariant solutions depending only on the variable y := x — Xt for 
some speed A and connecting two constant states (r_,M_) and (r+,M+) at infinity. The 
characterization of these "admissible" discontinuities is based on kinetic relations. (For 
background, see [20] and the references cited therein.) 

For instance, one important model of interest in fluid dynamics describes liquid- vapor 
flows governed by van der Waals's equation of state. For pioneering mathematical works 
on van der Waals fluids we refer to Slemrod et al. [29,30,13], who investigated self-similar 
approximations to the Riemann problem. The concept of a kinetic relation associated 
with (undercompressive) nonclassical shocks or phase boundaries was introduced by Abe- 
yaratne and Knowles [1,2], Truskinovsky [31,32], and first analyzed mathematically by 
LeFloch [19]. Kinetic relations and nonclassical shocks were later studied extensively by 
Shearer et al. [18,28] (phase transition), LeFloch et al. [14,15,16,21,5] (traveling waves, 
Riemann problem, general hyperbolic systems), and by Bertozzi, Shearer, et al. [7,6] 
(thin film model), as well as Colombo, Corli, Fan, and others (see [10,11,12,13] and the 
references cited therein). 

1.2. Approximation of undercompressive waves 

The present paper built on existing numerical work done by the first author and his 
collaborators [14,15,22,8,9] and devoted to the numerical investigation of non-classical 
shocks and phase boundaries generated by diffusive and dispersive terms kept in balance. 
These papers cover scalar conservation laws, and systems of two or three conservation 
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laws arising in fluid dynamics (Euler equations) and material science. Entropy stable 
schemes were constructed, and kinetic relations were computed numerically. Numerical 
kinetic functions turn out to be very useful to to evaluate the accuracy and efficiency of 
the schemes under consideration. 

Numerical experiments were also performed on undcrcompressive shocks for the models 
of thin films [6,24,25] and phase dynamics [26]. These authors did not investigate the role 
of the kinetic relation for these models, and one of our aims in the present paper is to 
tackle this issue. 

Recall that one effective numerical strategy to compute under-compressive waves is 
provided by the Glimm and front tracking schemes, for which both theoretical and nu- 
merical results are now available [20,9]. Both schemes converge to the correct non-classical 
solutions to any given Cauchy problem. The main feature of these schemes is to com- 
pletely avoid spurious numerical dissipation or dispersion, which this guarantees their 
convergence to the physically meaningful solution selected by the kinetic relation. How- 
ever, this class of schemes has some drawbacks: they are limited to first-order accuracy 
and require the precise knowledge of the Riemann solver, while numerical solutions may 
exhibit a noisy behavior. 

In contrast, in the last twenty years, modern, high-order accurate, finite difference 
techniques of (classical) shock capturing have been developed to deal with discontinuous 
solutions of hyperbolic problems. Extending these techniques to compute nonclassical 
shocks turned out to be quite difficult, however. Our aim in the present paper is to pur- 
sue this investigation of the interplay between dissipative and dispersive mechanisms in 
hyperbolic models, and to better understand how they compare with similar mechanisms 
taking place in finite difference schemes. 

This problem was first tackled by Hayes and LeFloch [14,15] where the importance of 
the equivalent equation associated with a scheme was pointed out and kinetic functions 
for schemes were computed. Next, the role of high-order, entropy conservative schemes 
was emphasized and further kinetic functions were determined numerically [22]. In the 
situation just described, the numerical solutions usually contain mild oscillations, which 
vanish out when the mesh is refined. This feature is entirely consistent with the behavior 
of traveling wave solutions to the underlying dispersive model. In consequence, total vari- 
ation diminishing (TVD) techniques should not be used to compute nonclassical shocks. 
Observe also that the diffusive-dispersive model serves only to provide a mechanism to 
select "admissible" solutions of associated hyperbolic eciuations. It has been established, 
for many hyperbolic models, that physically meaningful solutions can be characterized 
uniquely by pointwise conditions on shocks (Rankine-Hugoniot jump conditions, an en- 
tropy inequality, a kinetic relation, and, possibly, a nucleation criterion). 

Finite difference schemes can also be studied for their own sake, and one important 
issue is whether criteria can be found for a given scheme to generate nonclassical shocks 
or not. Some heuristics were put forward in [15] and allow one to distinguish between 
the following cases: 

(1) Schemes satisfying (a discrete version of) all of the entropy inequalities. In the context 
of scalar conservation laws, this is the case of the schemes with monotone numerical flux 
functions. For systems of two conservation laws (which admit a large family of (con- 
vex) mathematical entropies) this class includes the Godunov and the Lax-Friedrichs 
schemes (for small data, at least). These schemes, if convergent, must converge to a 
weak solution satisfying all entropy inequalities which, consequently, coincides with the 
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classical solution selected by the Oleinik or Kruzkov entropy conditions (scalar equa- 
tions) and the Wendroff or Liu entropy conditions (systems of two or more conservation 
laws) . 

(2) Schemes satisfying a single (discrete) entropy inequality but apphed to a system with 
genuinely nonlinear characteristic fields. The classification of Riemann solutions given 
in [16] shows that, again, only the scheme can converge to the classical entropy solution, 
only. 

(3) Schemes satisfying a single (discrete) entropy inequality but applied to a system with 
non-genuinely nonlinear characteristic fields. To decide whether such schemes are ex- 
pected to generate nonclassical shocks, one should determine the equivalent equation. 
One may truncate the equivalent equation and keep only the first two terms. An analy- 
sis of the properties of traveling waves for this continuous model provides an indication 
of the expected behavior of the scheme. Interestingly, the behavior depends on the sign 
of the dispersion coefficient and the sign of the third-order derivative of tlw^ flux. For 
instance, for first-order schemes applied to scalar equations one typically obtains, after 
further linearization in the neighborhood of the origin in the phase space, 

Nonclassical shocks have been observed when the flux is concave-convex (that is, v^) 
and a is positive, or else when the flux is convex-concave (that is, —v'^) and the 
coefficient a is negative. 



1.3. Purpose of this paper 

We will demonstrate here that the existence of under- compressive shocks and several 
typical behaviors of these nonlinear waves, especially the existence of a kinetic function, 
axe properties shared by many examples arising in continuum physics. To make our point, 
we present a number of physical models describing nonlinear wave dynamics dynamics: 
cubic flux, thin liquid films, generalized Camassa-Holm, van der Waals fluids. Specif- 
ically, we prove that kinetic functions can be associated to each of these models and 
we study their monotonicity, dependence upon (viscosity, capillarity, mesh) parameters, 
and behavior in the large. We uncover several new features of the kinetic function that 
have not been observed theoretically via analytical methods yet. It is our hope that the 
conclusions reached here nunicric;ally will motivate further theoretical developments in 
the mathematical theory of non-classical shocks. The work also provides further ground 
that not a single theory of entropy solutions but, rather, many theories of shock waves 
are required to accurately describe singular limits of hyperbolic equations, as supported 
by the framework developed in [20,21,23]. 

The outline of the paper is as follows. In Section 2, we present the physical models 
of interest and discuss briefly their analytical properties. In Section 3, after introducing 
some background on non-classical shocks and kinetic relations, we investigate the role 
of the equivalent equation. In Section 4, we estabhsh the existence of kinetic functions 
associated with each of the models and investigate their properties. In Section 5, we 
investigate van der waals fluids. Finally, Section 6 contains concluding remarks. 
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2. Models of interest 



We begin with a brief presentation of a few nonlinear hyperbolic models arising in 
continuum physics. 

2.1. Cubic conservation law 

It will be convenient to start with an academic example consisting of a conservation 

law whose flux- function admits a non-dcgcncratc inflection point. For simplicity and with 
little loss of generality as far as the local behavior near the inflection point is concerned, 
we can assume that the flux is a cubic function. After normalization, we arrive at the 
cubic conservation law 

ut + {u^)t = 0, u = u{t, x) e lR,t>0. (2.1) 

Wc arc interested in (discontinuous) solutions that can be realized as limits of diffusive- 
dispersive solutions of 

ut + {u^)x = (^Uxx + u^xx, u = ul^{t,x), (2.2) 

where a is a fixed parameter and e ^ 0. Recall that shock waves of (2.1) are solutions 
containing a single propagating discontinuity connecting two states U- , u+ at the speed 
A. These constants must satisfy the Rankine-Hugoniot relation 

A = — = u_ + U- u+ + u,, 

U+ — U- 

as well as the entropy inequality associated with the quadratic entropy U{u) := 

-X{ul-u^_) + ^{uX-ut)<0. 

As is now well-known, the solutions Ua ■= lime^o^^ may contain both classical shock 
waves satisfying the standard compressibility condition (equivalent here to the entropy 
criteria introduced by Lax) 

3 = /'(u_) > A > f{u+) = 34, (2.3) 

as well as non-classical shock waves which turn out to be under-compressive 

A < f{u±). (2.4) 

The characterization of the solutions generated by (2.2) as e ^ is provided as follows. 
Based on an analysis of all possible traveling wave solutions of (2.2), one can see that, 
for a given viscosity/capillarity ratio a and for every left-hand state U-, there exists a 
single right-hand state 

u+ = ^i{u_), (2.5) 

that can be attained by a non-classical shock. The function ip'^ is called the kinetic func- 
tion associated with the model (2.2). The existence of the kinetic function has been estab- 
lished theoretically for a large class of flux-functions and nonlinear diffusion-dispersion 
operators, including (2.2). The results are often stated in terms of the shock set, Sq(u^), 
consisting of all right-hand states u+ that can be attained from a given left-hand state 
U- by a classical or by a non-classical shock. Note also that instead of the relation (2.5), 
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one can equivalently prescribe the entropy dissipation of a non-classical shock, that is 
the kinetic relation can be expressed in the following form observed in [19]: 

- A {U{u+) - U{u-)) + F{u+) - F{u-) 

= - [ U"{v{y))vy{yfdy<0, ^^"^^ 

where y i— > v{y) denotes the traveling wave trajectory connecting u- to u+. In orther 
words, the entropy dissipation must be prescribed on a non-classical shock. As far as 
the specific model with cubic flux and linear diffusion and dispersion is concerned, the 
kinetic function and the shock set can be expressed explicitly by analytical formulas as 
follows. The kinetic function associated with (2.2) reads 



-u- — 5/2, u- < —a, 

-U-/2, \u^\<a, (2.7) 

-U- + a/2, u- > a. 



with a := ^/{8/3a), while the corresponding shock set is 

{{u-,a/2]u{-U--a/2}, U- < -a, 
[—u-/2,u-), —a<U-<a, 
{-u^ +a/2}u[-a/2,u^), U- > a. 

Observe that (p^ converges to —u/2 when a ^ 0, and that the shock set converges to 
the standard interval [—u/2, u] determined by the Oleinik entropy inequalities. See [20] 
for details. 

2.2. Thin liquid films 

The thin liquid film model (e, a being positive, scaling parameters) 

ut + {u^ - u^)x = e {u^ Ux)x - a(? {u^ Uxxx)x (2.8) 

describes the dynamics of a thin film with height u = u{t, x) moving on an inclined flat 
solid surface. The (non-convex) flux-function 

f{u) = u^ -u^, < u < 1, 

represents the competing effects of the gravity and a surface stress known as the Marangoni 
stress. The latter arises in experiments due to an imposed thermal gradient along the 
solid surface. The fourth-order diffusion is due to surface tension, and the second-order 
diffusion represents a contribution of the gravity to the pressure. This model was in- 
troduced and extensively studied by Bertozzi, Miinch, and Shearer [6]. The existence of 
non-classical traveling waves was established analytically in [7], and various numerical 
studies were performed [25,24] which exhibited subtil properties of stability and insta- 
bility of these waves. For more analytical background, see also the convergence theory 
developed by Otto and Westdickenberg [27]. 

Prom the standpoint of the general well-posedness theory the kinetic relation is an 
important object which is necessary to uniquely characterized the physically mcaiiiiigful 
solutions. A kinetic function was not exhibited in the work [7] which, instead, relied 
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on non-constructive arguments. Prom the existing literature, it is not clear whether a 
concept of a kinetic relation could be associated with the equation (2.8) and, if so, and 
whether such a kinetic function would enjoy the same monotonicity properties, as the 
ones established earlier for conservation laws regularized by viscosity and capillarity. This 
issue will be addressed in Section 4. 

Observe that the flux f{u) here is positive on the interval (0, 1), increasing on (0, 2/3) 
and decreasing on (2/3, 1). It admits a single inflection point at u = 1/3. To every point 
u e (0, 1/3) we can associate the "tangent point" (p{u) G (1/3, 1/2) characterized by 

/'(.■(»)) = M^4fM, 

u — (p»{u) 

or 

V>\u) := . 

The same formula maps also the interval (1/3, 1) onto (0, 1/3). This function allows us 
to dcflnc classical shock waves associated with the equation. When u G (0, 1/3) the left- 
hand state u can be connected to the right-hand state ip^ {u) by a contact discontinuity. 
When u € (1/3, 1), (p\u) is the left-hand state and u is the right-hand state. 
Another important function is provided by considering the entropy dissipation 

D{u.,u+) := -^{ul - ) + ^(4 - ut) + liuX - ut). 
where the shock speed A is given by 

A := - =u+ + u.- {ul + u+u. + ut.). 

U+ — U- ^ 

The zero dissipation function (p^ is by definition the non-trivial root of D, i.e. 

Diu,ipi{u))=0, 

or 

V"(w) •= 3 

The function (/j" maps the interval (0,2/3) onto itself. 

According to the theory in [20] the range of the nonclassical shocks is limited by the 
functions (f^ and (p". Precisely, for a nonclassical shock connecting it_ < 1/3 to u+ > 1/3, 
the right-hand state must have 

The sign are reversed for a decreasing nonclassical shock. 
2.3. Generalized Camassa-Holm model 

We consider a generalized version of the Camassa-Hohn equation 
Ut + f{u) 



(2.9) 



7 



which arises as an asymptotic higher-order model of wave dynamics in shallow water. 
The second-order and third-order terms are related to the viscosity and the capillarity 
of the fluid. When the flux / is nonconvex, for instance 

f{u) = 

nonclassical shocks may in principle arise. We will demonstrate in this paper that indeed 
solutions may exhibit nonclassical shocks and establish the existence of an associated 
kinetic function. 



2.4. Van der Waals fluids 



Compressible fluids are governed by the following two conservation laws: 

dtT — dxU — 0, 

^ ^ ^ 2. (2-10) 

Here, u and r represent the velocity and the specific volume of the fiuid, respectively, while 
a is a non-negative parameter representing the strength of the viscosity. The pressure 
law p = p{t) is a positive function defined for all r € (0, +oo) and of the following van 
der Waals type: there exist < a < c such that 

p"{t)>0, r e (0,a) U (c,+oo), 

p"{t)<0, re(a,c), (2.11) 
p'{a) > 0, 

and 

lim p{t) = +O0, lim p{t) = 0. (2.12) 

T — >Q T — * + 00 

The left-hand side of (2.10) forms a first-order system of partial differential equations, 
which is of elliptic type when r belongs in the interval {d, e) characterized by the con- 
ditions 0<d<a<e<c and p'{d) = p'{e) = 0. It is of hyperbolic type when 
T G (0,d) U (e,-|-oo) and admits the two (distinct, real) wave speeds ±-\/— p'(t). 



3. A conjecture on the equivalent equation 

3.1. Kinetic functions associated with difference schemes 

We now begin the discussion of the numerical approximation of the solution Ua gen- 
erated by (2.2). The discussion applies to general conservation laws of the form 

dtU + dxf{u) = 0, 

where the fiux-function / admits a single infiection point. To any finite difference scheme 
associated with (2.1) one can in principle associate a kinetic function, which we denote 
by V'a- It may seem natural to request that Va coincides with the kinetic function tp^ 
associated with the given model. However, it has been observed by Hayes and LeFloch 
[15] that, at least for all finite difference schemes that have been considered so far, 

(3.1) 
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This discrepancy is due to the fact that the dynamics of non-classical shocks is determined 
by small-scale features of the continuous model (2.2) which can never be fully mimicked 

by a discrete model. 

Given this perspective, the next natural question is to compare the kinetic functions 
tp^^ and Wc require that a scheme be a good approximation of the given, continuous 
model in the following sense. Its equivalent equation obtained by formal Taylor expansion 
should have the "correct" form 

Vt + {v^)x = cVa;a; + ae'^Uxxx + 0{h''), h = ce, (3.2) 

where c is a constant and q represent the order of accuracy of the scheme. Here, v denotes 

the numerical solution and h denotes the discretization parameter. It is important to 
observe that all of the schemes considered in the present paper are first-order accurate, 
only, as far as the hyperbolic equation (2.1) is concerned. They are, however, high-order 
approximations of the augmented model (2.2). For clarity, wc emphasize the dependence 
in q and denote by ^/i^'^ the kinetic function associated with a scheme whose equivalent 
equation is (3.2). 

L(rt us introduce in this paper the following: 
Conjecture 3.1 As q oo the kinetic function V'a^ associated with a scheme with 
equivalent equation (2.7) converges to the exact kinetic function v^^j 

lim VL'' = (3.3) 

q — ^oo 



Rigorous results pointing toward the validity of this conjecture can be found in [4] 
which studies the role of relaxation terms in traveling wave solutions of conservation 
laws. A closely related problem was tackled by Hou and LeFloch [17] who considered 
nonconservativc sc;hemcs for the computation of nonlinear hyperbolic problems. In these 
problems, small-scale features are critical in selecting shock waves, and the equivalent 
equation have been found to provide a guide to designing difference schemes. We will 
not try here to support the above conjec;ture on theoretical grounds, but we propose 
to investigate it numerically. As we will see, very careful experiments axe necessary. 
We consider a large class of schemes based on standard differences, and obtained by 
approximating the spatial derivatives arising in (2.2) -with e replaced by the mesh size 
h (up to a constant multiplicative factor) 



f{u)a;, hUxx, h 



2 , 



by high-order finite differences, so that the overall scheme is of order q at least. For 
completeness we list below the corresponding expressions. We denote by Xj the points of 
spacial mesh and by Ui the approximation of the solution at the point Xi. We also use 

the notation fi := f{ui). Based on the above we arrive at semi-discrete schemes for the 
functions Ui = Ui(t). For instance, using fourth order discretizations above we obtain 



-If 

e / 1 4 5 4 1 \ 



h 
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Observe that this is a fully conservative scheme, in the sense that is indepen- 

dent of t (assuming, for instance, periodic boundary conditions). To actually imple- 
ment the above algorithm, we use a Runge-Kutta scheme. Defining U{t) = {ui{t)) (with 
i= ...,—1,0,1,...), the semi discrete scheme takes the form 

^{t) = R[U{t)], (3.5) 

where R[U{t)] represents the spatial discretization. This system of ordinary differential 
equations is solved numerically by employing an s-stage Runge-Kutta scheme defined as 
follows 

fc-i 

g>' = R(u^ + AtJ2ak,j9'), 

(3.6) 

Ijn+l ._ + Ai ^fefe/. 

fc=l 

For example, the non-zero coefficients of the fourth-order Runge-Kutta scheme are given 
by 

^2,1 = 1/2, 03,2 = 1/2, 04,3 = 1, 

&i = 1/6, 62 = 1/3, &3 = 1/3, 64 = 1/6. ^ ■ ' 

Coefficients of a sixth and an eighth-order Runge-Kutta scheme are found in [33] and 
[34], respectively. 



3.2. Numerical experiments 



We now determine numerically the kinetic function associated with the schemes de- 
scribed in the previous section. For each q = 4,6,8,10 and for selected values of the 
parameter a we compute the function tp^'^. For the cubic flux function, two typical types 
of non-classical waves arise, which are under-compressive shock followed by a rarefaction 
wave (Figure 1, left), and a double shoc;k struc;ture (Figure 1, right). In order to plot a 
single kinetic function it is necessary to solve a large number of Riemann problems for 
various values of left-hand state U- and to ensure that the right-hand state is picked up 
in an interval where a non-classical shock does exist. In the numerical experiments it is 
often convenient to select the right-hand state so that the Riemann solution contain two 
shocks. According to the construction algorithm in [21] the middle state between the two 
shocks is therefore the kinetic value = V'^ ''(«-) • 

Several difliculties arise. First, when u_ is close to the origin all the waves become very 
weak and it is numerically difficult to identify the middle state. Second, in a range of 
parameter values, the two waves may propagate with speeds that are very close and this 
again makes difficult the computation of the middle state. Finally, solutions do contain 
mild oscillations, especially for small values of a and this again introduces some numerical 
error. Due to these constraints we need to use a rather fine mesh, with h of the order of 
1/1000. 

The following plots allow us to investigate numerically the convergence of the kinetic 
function toward the exact kinetic function, (2.7), which is a piecewise affine function 
of the variable U- . 
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Fig. 1. Typical cases of non-classical waves for the cubic flux function. 

In the following numerical tests, the CFL number was taken to be as large as possible 

in each run and was identical for all schemes (4th order to 10th order). It was observed 
that increasing the order of Runge-Kutta scheme (in the temporal integration) from four 
to six and eighth practically does not change the results, therefore, a tenth order Runge- 
Kutta scheme was not examined. On the other hand, the order of spatial discretization 
was found to be very important and effectively made an important difference. 

First of all. Figure 2 shows the right-hand state u+ versus the parameter a for different 
schemes. Here, we have used h = 0.001, u_ = 10, and e = 5/i. As a increases, the results 
of different order numerical schemes become closer and closer to the exact solution. For 
large values of a, e.g. a > 10, even the fourth-order scheme gives satisfactory results. 
Since the solution of (2.1) is, in fact, the limit of diffusive-dispersive solutions of (2.2), 
with a fixed and e ^ 0, we conclude that (provided a is sufficiently large) even the 
fourth order scheme converges to the exact solution. Figure 2 shows also that, for a fixed 
value of a (sufficiently large), as the order of accuracy increases, the numerical solution 
converges to the exact one. These results strongly support Conjecture 3.1. 

Next, Figure 3 (left column) shows u+ versus U- for a=l, 4 and 6 respectively, with 
h = 0.005, and e = 5h. A grid of 2000 nodes was used in all runs. Again, by increasing 
the order of accuracy, the numerical solution converges to the exact one, for sufficiently 
small left-hand state m_. However, the suitable value of a depends on u_ and increases 
as a does. It should be mentioned that very high values of a (which are in need for 
large U-), are not satisfactory since they lead to high oscillatory results (recall that a is 
dispersion-diffusion ratio). It was also observed that for small values of a, high frequency 
oscillations occur before the non-classical shock, while for large values of a, low frequency 
oscillations take place after the rarefaction (not shown). 

Finally, Figure 3 (right column) shows the scaled entropy dissipation (j){s)/s'^ versus 
the shock speed s for a = 1, 4, 6, respectively. We use the quadratic entropy U{u) = v? jl 
in computing the entropy dissipation. In terms of u_ and u+, this is 

(^(■u_, u+) = (u+ — u_)^(m^ — u?_), (3.8) 
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Fig. 2. uj^ versus a for different schemes. 

and from the Runkine-Hugoniot relation, the shock speed is 

s — + u^u+ + u\. (3-9) 

The same feature is observed as before and by increasing the order of accuracy; the 
numerical entropy dissipation converges to the exact one (provided a is sufficiently large) , 
which is in accordance with Conjecture 3.1. 



3.3. Spatial discretization parameter versus diffusive-dispersive parameters 

To conclude this section let us set h = ce and plot the corresponding numerical results 
with the sixth-order scheme and some fixed values of h and e. It was observed that when 
c is either too small or too large, the numerical results deteriorate and an approximate 
kinetic function can no longer be associated with the numerical approximations. 

As shown in Figure 4 for typical values m_ = 10, a = 1 and CFL=0.5, the specific 
choice of c does not change the general feature of the results. However, as we increase 
the value of c, all schemes converge to the exact solution and the tenth order scheme 
converges faster. Moreover, it is concluded from Figure 4 that the sixth order scheme 
with c = 10 is a reasonable choice in terms of accuracy and computational efficiency. For 
this value of c (that is, e — lO/i) with h — 0.005, Figure 5 shows the kinetic function 
Uj^ versus u„ (left) and the corresponding scaled entropy dissipation (f>{s)/s'^ versus the 
shock speed s for a = 1,4,6 (right) which are in accordance with the previous results 
about the convergence of the numerical results to the exact one. 
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M+ versus U- 



(j){s)/s^ versus s 
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4. Fourth-order models 

We now return to the examples of Section 2 and discuss the thin film and Camassa- 
Holm models. We will show that solutions of the Riemann problem may admit under- 
comprcssivc shocks, and wc will determine the corresponding kinetic function in a sig- 
nificant range of data. Our results in the previous section indicate that we can continue 
this investigation with six-order schemes, which we will do in (most of) the forthcoming 
tests. 

4.1. Thin liquid films 

Consider the model of thin liquid films together with the physical rcgularization terms. 
Our goal is to (numerically) compute a kinetic function associated with this model. It 
should be reminded that the physically relevant range for the solutions is the interval 
(0, 1). This model has a flux with a single inflection point, as is the case with the cubic 
conservation laws, and the interesting aspect is the particular form of the rcgularization. 
In view of experiments performed in [7], it may be anticipated that a unique Riemann 
solver can not be associated to a given set of initial data, and that the specific initializa- 
tion adopted in implementing a given scheme may play a role. We will first show that a 
kinetic relation can be associated to a family of initial data and a given approximation 
of these initial data. 

Since the fiux has a convex-concave shape rather than a concave-convex form (contrary 
to the cubic fiux function) , the nonclassical shock is fast undcrcompressive and it is more 
convenient to draw the kinetic functions from right to left, that is. 

The calculations here are more delicate than those performed with the third-order rcgu- 
larization, since we are now dealing with a fourth-order rcgularization. Therefore, for the 
discretization of the terms arising in (2.8), a scheme must be fifth order in accuracy at 
least. We thus exclude the fourth-order scheme used earlier and we concentrate attention 
on the sixth-order scheme. 

It should be mentioned that the (fourth-order) rcgularization is singular at u = 0, so 
tests involving values close to u = are challenging and the kinetic function should have 
some degenerate behavior as some values in the solution approache u = 0. 

The thin film model (2.8) may be written in the conservative form 

Wt+S'(«)x = 0, (4.1) 

with 

g{u) := -u^ - {S - u^xx), (4.2) 
where we have set e = S and ae^ = 1. 

The numerical solution of (4.1) is obtained by calculating the terms and u^xx in 
(4.2) using the sixth order scheme (Appendix A), computing the nodal values of g, and 
finally calculating Qx again using the sixth order scheme (Appendix A). The numerical 
experiments are performed here with S = 0.1 /i and h = 1 using a computational grid of 
1000 nodes and the initial condition used in [7] 

u{x) = (tanh(-x + 100) + 1) ~ + ur. (4.3) 
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Fig. 6. Wave structure for uji = 0.1 and u^, = 0.5 and 0.6. 

The numerical solutions of (4.1) for ur = 0.1 and ul = 0.5 and 0.6 with At = 0.6481 
at time t = 1037 are shown in Figure 6. For the case ul = 0.5 a double shock structure 
is observed while for ul = 0.6 a rarefaction wave and an under- compressive shock are 
generated. The speed of the non-classical shock is similar for the two cases, as expected. 

In agreement with what was observed (and proven rigorously in [5]) for the diffusive- 
dispersive regularizations, when 0.2 < ur < 0.4, i.e. around the inflection point u= 1/3, 
only classical waves are observed as shown e.g. in Figure 7 for ur = 0.3 and different 
values of Ml (with At = 0.388 at time t = 698). For 0.4 < ur non-classical waves arise 
even when ul is small and bcc;omc c;k;arly visible as ur increases as shown in Figure 8. 

Figure 9 shows the numerical solutions obtained for At = 0.37, at time t = 2148 for 
Ur = 0.6 and two cases ul = 0.1 and ul = 0.2. For the case ul = 0.2 a double shock 
structure is observed while for ul = 0.1 a rarefac;tion wave and an under-comprcssivc 
shock are generated. It should be mentioned that when Ul is selected large enough, only 
classical waves are observed as shown for ul = 0.5 and 0.68 (with the same data) in 
Figure 10. 

Varying the left-hand and right-hand initial states and keeping fixed all other param- 
eters, for 6 = T]h with r] = 0.1 we obtain the numerical kinetic function (u_ as a function 
of Ur) shown in Figure 11 (left). Interestingly enough, this function is decreasing, as 
required in the general theory of nonclassical shocks. It is also observed that the nu- 
merical kinetic functions converge as the order of the accuracy of the numerical scheme 
is increased. Note that, for sufficiently large riglit-liand values of ur, negative left-hand 
values U- are obtained which are physically irrelevant. The scaled entropy dissipation 
(t){s)/s'^ versus the shock speed s is also shown in Figure 11 (right). 

Finally, we study the effect of ry = 5/h for Uj. — 0.8. Figure 12 presents a sample of the 
solution with various schemes for r] = 4.7. In Figure 13, it is shown that by increasing 
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Fig. 7. Wave structure for un = 0.3 and various values of uj^. 




Fig. 8. Wave structure for uj^ = 0.05 and ujj = 0.4, 0.5 and 0.6. 

7], the numerical results converge to the limiting value with the tenth order scheme. The 
optimum value of 77 depends on Ur and it is increased as Ur does (not shown) . In Figure 
14 we have plotted the kinetic functions for various schemes where a variable (optimum) 
T] (depending on Ur) is employed. The convergence predicted in Conjecture 3.1 is clearly 
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Fig. 9. Wave structure for un = 0.6 and u/, = 0.1 and 0.2. 
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Fig. 10. Wave structure for un = 0.6 and u/, = 0.5 and 0.68. 



confirmed here. 

In order to investigate the effect of the initial conditions, we consider the following 
initial condition 
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Fig. 11. Kinetic function (left) and scaled entropy dissipation (right), for rj = 0.1 based on the initial 
conditions (4.3). 
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Fig. 12. A sample of the solution with various schemes for Ur = 0.8 and rj = 4.7. 

' (0.35 -UL)tanh(a;- 80)/2 + (0.35 + ml)/2, x< 105, 



{uR - 0.35) tanh(a; - 130)/2 + (0.35 + ur) /2, x > 105. 



(4.4) 



The numerical solutions of (4.1) for ur^ = 0.7 and = 0.05 with At = 0.27 and Ax = 1 
at time t = 18792 are shown in Figure 15. As shown in Figure 15, the two choices of 
numerical solutions lead to different results. However, such a dependence on the initial 
data is only limited to few cases. 
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Fig. 13. U— versus r] = 5/h for Ur = 0.8. 
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Fig. 14. Kinetic function for thin films, based on the initial conditions (4.3) with optimum rj. 

4.2. Generalized Camassa-Holm model 

It is expected that the kinetic function of the Camassa-Holm model is qualitatively 
similar to that of the diffusive-dispersive model. However, some differences, at least quan- 
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Fig. 15. Numerical results obtained using the initial conditions (4.3) (left) and (4.4) (right). 

titatively are expected. In the following, we draw a kinetic function for the Camassa-Holm 
model and compare with the diffusive-dispersive model with the same coefficients. First, 
we investigate the behavior of the function for large values of u. 

Figure 16 shows the results obtained with h = 1, e — O.OOBh, a = 1 and CFL=0.5 
using a computational grid of 1000 nodes. The kinetic function for the Camassa-Holm 
model (Figure 16, right) is larger than the linear diffusion-dispersion case (Figure 16, 
left). 

The comparison between the Camassa-Holm model and the linear diffusion-dispersion 
model allows to see that the kinetic function of the two models is similar for small values 
of u but different for large values of u. Typically, for large values of u, the Camassa- 
Holm model leads to larger kinetic values. This behavior is due to the nonlinearity of the 
regularization in the right-hand side of (2.9). It would be interesting (and challenging) 
to apply the techniques in [5] and prove the existence and monotonicity of the kinetic 
function associated with the Camassa-Holm model. 

We then consider small values of u. Figure 17 shows the results obtained for both 
models with h = 1, e = O.lh, a = 4 and CFL=0.5 using a computational grid of 1000 
nodes. Contrary to the case of large values of u, the behavior of the function for small 
values of u is the same for both models. 



5. Kinetic functions associated with van der Waal fluids 

Our next objective is to investigate the properties of non-classical solutions to hyper- 
bolic conservation laws whose flux-function admits two inflection points. We consider 
the phase transition model (2.10) presented in Section 2. As we will see, the dynamics 
of non-classical shock waves for general flux is much more intricate than the now well- 
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Fig. 17. The kinetic function (small values of u) for the linear diffusion-dispersion case and Camassa-Holm 
model. 



understood case of a single inflection point. In particular, we demonstrate numerically 
that the Riemann problem may well admit arbitrarily large number of solutions. This 
happens only when the flux admits two inflection points at least and provided the disper- 
sion parameter is sufficiently large. We then determine the kinetic function and observe 
its lack of monotonicity. 

We consider an equation of state having the same shape as that of van der Waals 
fluids, described by the (normalized) equation 

1 3 

P(^) ■= (3r-l)i+i/C "7^' > ^^-^^ 
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for some positive constant C = 1/(7 — 1) where 7 G (1,2). 

To exhibit the dynamics of non-classical shocks for general flux, we solve the Riemann 
problem for left- and right-data; 

(i) First, we determine the wave structure of the solution for each fixed left-hand state, 
that is, we identify the waves (classical/non-classical shock or rarefaction) within 
the Riemann solution. In turn, by varying the left-hand state, we identify regions in 
the plane in which the structure of the Riemann solution remain unchanged as we 
change the Riemann data. This provides us with a representation of the Riemann 
solver in the plane. 

(ii) Second, within the range of Ui where non-classical shocks leaving from ui arc avail- 
able, we determine the corresponding kinetic functions that are needed to char- 
acterize the dynamics of non- classical shocks. It is expected that more than one 
kinetic function will be needed in some range of ui, but a single kinetic function 
maybe sufficient in certain intervals. 

To begin with, we consider a simplified case 

dtT — dxU = 0, 

dtu + dxP{T) = e dxxU, 

with 

/ ^ 3 /r ON 

and R = ^ and T = 1.005. This flux function has two inflection points at r = 1.00996 
and 1.8515 and it is shown in Figure 18. 

In order to investigate possible wave cases, various test cases are performed by keeping 
three of four variables {ul, tl, tr and u^) constant and changing the fourth one. Some 
typical wave structures for ul = 0.35, — 0.8, = 2 and ujj = 0.5,1.5 and 2 are 
shown in Figures 19, 20 and 21. 

Next, we select ur = 1,tl = 0.8, tr = 2 and ul variable from to 2. Note that the two 
inflection points are between the selected tl and tr. We study the kinetic function for 
this case. Numerical results are obtained using a c;oniputational grid of 2000 nodes and 
h = 0.0005 with e = 0.00003. The wave structures for this case are given in Figures 22 
and 23 (corresponding to the case ur = 1). 

The kinetic function obtained for ur = 1 is shown in Figure 24, which shows the 
relation between right and left states of the nonclassical shock for r. Note that the 
kinetic functions are obtained only for ul < 1-2, where clear nonclassical shocks are 
generated. It is observed that for large w_, all schemes give close results, while for small 
the results largely depend on order of accuracy of the numerical scheme. Again, as the 
order of accuracy increases, the numerical solutions appear to converge, which supports 
Conjecture 3.1. The kinetic function in this case is monotonic; and strictly decreasing. 

In order to check whether the kinetic function is single-valued, we now consider ur = 
1.5. The kinetic function obtained in this case for r is shown in Figures 25. As it it 
observed, the results axe similar to those of kinetic function in Figure 24. Similarly, the 
cases with UR = Q,UR = i and ur = ^ were also tested, and they led to the same kinetic 
function. This shows that the kinetic function is single- valued here, at least for the cases 
considered above. 

Next, we add the capillarity effects in (5.2) 
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Fig. 18. Equation of state considered in (5.3). 
dtT - dxU = 0, 



(5.4) 



We repeat the experiments with the same data (i.e. uji = I, tl = 0.8, th — 2 and varying 
Ml), but now including capillarity effect with a = 1. The capillarity coefficient is ae^ as 
in the previous cases. The kinetic functions for various schemes are shown in Figure 26, 
(which shows the relation between right and left states of the nonclassical shock for r) , 
and they are compared with the case of no capillarity effect. As it is observed, for a given 
T_, the capillarity effect leads to smaller values of r+. 

5.1. A piecewise linear pressure function 



In order to highlight the effect of a pressure function with two inflection points, here 
we use a piecewise linear flux function which was already studied in [3] , given by 



p{t) 



-7t+ 10, 
4r- 1, 

-^+12, 
T 4 

^r+15' 



T< 1, 

1 < T < 2, 

2 < T < 4, 

4 < T, 



(5.5) 



and shown in Figure 27. In order to explore various regimes occurring by this pressure 
function, we set — 0.9 and tr = A. The two inflection points are again between the 
selected tl and th. 
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Note that, since the pressure function is piecewise hnear, r in the system 5.4 does not 
depend on the specific choice of ul and un, but depends only on their difference uj^ — ur. 
This is easily checked by a change of variable u + c, where c is any constant. Hence, 
in the numerical experiments without loss of generality, we can fix Ufl = 1, and change 
Ul only. 

Numerical results are obtained using a computational grid of 2000 nodes and h = 
0.0005 with e = 0.001 and a = (no capillarity effects). 

Three regimes are identified, as explained now: 

- Regime A: lA < ul and a stationary shock. 

In this regime, a stationary shock wave is generated in the center. A typical solu- 
tion is plotted in Figure 28 for ul = 1.5. The kinetic function in this regime based 
on the fourth order scheme is shown in Figure 29, and it appears to be a linear function. 

- Regim,e B: — 1.2 < < 1-4 and a non- stationary shock. 

In this regime, a (left-going) nonclassical shock wave is generated. A typical picture 
of this regime is shown in Figure 30 for ul = 1- As decreases, the shock speed 
is increased, and the left-hand and right-hand values of the nonclassical shock reach 
some limiting values. The kinetic function in this regime is non-monotone and not 
single-valued. 

- Regime C: Ul < —1.2 and a non-stationary shock with fixed left- and right-hand values. 
In this regime, a (left-going) nonclassical shock wave is generated, but the values r_ 
and T+ are fixed to the limiting values of Regime B awhile the specific volume r is 
increased on the right-hand side of the nonclassical shock and generates a right-moving 
wave. A typical solution for this regime is plotted in Figure 32 for = —2. No kinetic 
function arises in this regime since t_ and t+ do not change. 

Finally, kinetic functions for both Regimes A and B are given in Figure 33, using 
fourth- and eighth-order schemes. Note that the kinetic functions axe almost identical 
in the Regime A (stationary shock) for both schemes, but they are largely distinct in 
the Regime B. 



6. Concluding remarks 

The specific contributions made in this paper are three-fold. 
- First, we have investigated the role of the equivalent equation associated with a scheme. 
We stated a conjecture and provided numerical evidences demonstrating its validity. 
We have shown that the kinetic function associated with a finite difference scheme 
approaches the (exact) kinetic function derived from a given (viscosity, capillarity) 
regularization. Interestingly, the accuracy improves as its equivalent equation coincides 
with the diffusive-dispersive model at a higher and higher order of approximation. We 
observed that the spatial accuracy plays a more crucial role than the temporal accuracy. 
In both the continuous model and the discrete scheme, small scale features are critical 
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Fig. 19. Typical wave structure for van dcr Waal fluids; left u and right r, for = 0.5. 
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Fig. 20. Typical wave structure for van dcr Waal fluids; left u and right r, for uj^ = 1.5. 



to the selection of shocks. Hence, the balance between diffusive and dispersive features 
determines which shocks are selected. These small scale features can not be quite the 
same at the continuous and at the discrete levels, since a continuous dynamical system 
of ordinary differential equations can not be exactly represented by a discrete dynamical 
system of finite difference equations. The effects of the regularization coefficients on 
the kinetic functions were also investigated. 

Second, we considered fourth-order models. We demonstrated that a kinetic function 
can be associated with the thin liquid film model. We also investigated a generalized 
Camassa-Holm model, and discovered nonclassical shocks for which we could determine 
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Fig. 21. Typical wave structure for van der Waal fluids; left u and rigfit t, for uji = 2. 
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Fig. 22. Wave structure (r), corresponding to the case = 0.8, th = 2, = 1 and variable m^. 

a kinetic relation. In both models, the kinetic function was found to be monotone 
decreasing, as required in the general theory [20]. 

Third, we investigated to what extent a kinetic function can be associated with van der 
Waals fluids, whose flux-function admits two inflection points. We established that the 
Riemann problem admits several solutions whose discontinuities have viscous-capillary 
profiles, and we exhibited non-monotone kinetic functions. 
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Fig. 23. Wave structure (u), corresponding to the case = 0.8, th = 2, un = 1 and variable uj^. 
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Fig. 24. The kinetic function for t obtained for un = 1. 

In summary, the new features of the kinetic functions observed here suggests directions 
for further challenging theoretical work on the well-posedness theory for the Cauchy 
problem associated with nonlinear hyperbolic problems. 
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Fig. 26. The kinetic function for t obtained for un = 1, and capillarity effect with a = 1, using the tenth 
order scheme. The dash curve shows the case without the capillarity effect. 
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Fig. 27. A piecewise linear flux function resembling the van der Waals pressure law. 
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Fig. 28. A typical wave structure in Regime A; u (left) and r (right) at time t = 0.12. 
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Appendix For completeness we list here the high-order discretizations used in this paper which, for 
instance, we state for a conservation law with conservative variable u and flux / = f{u). 

- Fourth order discretization (g = 4): 

,2 1 4 5 4 1 

h Uxx = - —Ui-2 + - -«i + - — "i+2, 

,3 1 1 

" Uxxx = —-Ui-2 + Ui-l — Ui+l + -Ui+2- 

- Sixth order discretization (q = 6): 

^^"^ ^ ~60"^'~^ 20'^*"^ ~ 4^^'^^ 4^^+^ ~ 20^^'^^ 60'^*'''^' 

1 3 3 49 3 3 1 

2 180 40 4 36 -1 + 40 + 180 + ' 

1 1 13 13 1 1 

y Uxx. = -Ui-3 - -«i_2 + - l + -U,+2 - -Ui+3. 
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4 

"i + -Ui+i 
5 



- Eighth order discretization (g = 8): 

1 4 14 4 1 

280 105 5 5 5 5 

105^'+' 280-^'+''' 
h? -I 4 1 4 205 

Uxx = Wi_4 H Wi_3 — — Ui—2 H — Ui 

2 1120 315 10 5 144 

14 1 

Ui+2 H Ui+3 Ui+4, 

10 ^ 315 ^ 1120 ^ 

/i^ -7 1 169 61 61 

Uxxx = + Uj — 'i — Uj — O + Uj — ] — U-i-\--\ 

6 '^'^'^ 1440 20 720 180 180 + 

169 1 7 

+ 720"^+^-20"'+^+14S0"^+^- 

- Tenth order discretization {q = 10): 

_-l 5 5 5 55 5 

"1260 ^ Wa^'^* ~ 84 + 21^*-' " 6^'-' ^ 6^'+' " 21^'+' 

5 5 1 

+ 84^'+' " 504-^'+* 1260^'+'' 
/i2 1 5 5 5 5 5269 5 

'<J'XX = fi— 5 Ui—n H Hi— 3 Ui—2 H Mi— 1 Ui H Uj+l 

2 6300 2016 252 42 6 3600 6 ^ 

5 5 5 1 

Ui+2 H Wi+3 «i+4 H Ui+5, 

42 ^ 252 ^ 2016 ^ 6300 ^ 
ftS 41 1261 541 4369 1669 

Uxxx — t*i — ~ Ui — 4 + W?— 3 ^ Uj — O + — \ 

6 36288 90720 6720 15120 4320 

1669 4369 541 1261 41 

— ^i-t-l H ^7^2 ~ ^?-t-3 H ^7+4 ^ ^7-t-5- 

4320 ^ 15120 ^ 6720 ^ 90720 ^ 36288 ^ 
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